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Abstract. We present a new multi-fluid, grid MHD code PIERNIK, 
which is based on the Relaxing TVD scheme ( |Jin and Xin, 1995] ) . The 
I . original scheme (see Trac & Pen (|2003|) and Pen et al. (|2003[) ') has been 

O ' extended by an addition of dynamically independent, but interacting 

fluids: dust and a diffusive cosmic ray gas, described within the fluid 
approximation, with an option to add other fluids in an easy way. 
The code has been equipped with shearing-box boundary conditions, 
and a selfgravity module, Ohmic resistivity module, as well as other 
facilities which are useful in astrophysical fluid-dynamical simulations. 
^ ' The code is parallelized by means of the MPI library. In this paper 

we present Ohmic resistivity extension of the original Relaxing TVD 
MHD scheme, and show examples of magnetic reconnection in cases of 
' uniform and current-dependent resistivity prescriptions. 

<N 

qq ■ 1 Resistive MHD equations 

o : 

Picrnik MHD code is capable of dealing with resistive MHD equations 
( |Pawlaszek, 2008) , which involve the resistive dissipation term in the induction 
■ equation 

a t B = Vx(vxB)-Vx(,J), (1.1) 

where J is the electric current density. Due to the conservative nature of the 
basic MHD algorithm, any decrement of magnetic energy during the update of 
magnetic field results in a corresponding increment of thermal energy. Therefore, 
no additional source term related to resistivity is needed in the energy equation. 

The resistivity algorithm of PIERNIK code relies on additional terms of re- 
sistive origin in the electromotive force, supplemented to the original Constraint 
Transport (CT) scheme implemented by Pen et al. (12003R within the Relaxing 
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TVD scheme. In the present version of PIERNIK we have implemented a current- 
dependent resistivity, according to Ugai l|1992p 

v = m+ m(J 2 - 4it)0(|J| - IJcritl), (1.2) 

where rji is a constant corresponding to uniform resistivity, 772 represents the 
anomalous resistivity, which switches on when the current density exceeds a certain 
critical value J cr itj and O is the Heaviside step function. 



2 Resistivity algorithm 

The relation between electric currents and magnetic fields in magnetized media is 
given by Ampere's law 

J = VxB, (2.1) 

where the code normalization of current densities is applied, to get rid of the 
factor c/Att in Ampere's law. To solve the Ampere's equation numerically in a 
grid scheme we set the location of current components at the cell edges, while 
the current-dependent resistivity values are initially placed at the cell corners 
(see Fig. [T]for details). To improve the stability of the resistive part of the MHD 




Fig. 1. Position of the magnetic field B components, the electric current J components 
and the magnetic resistivity 77 on a computational cell. Braces indicate interpolation of 
77 to the position of electric current density J 

scheme a weighted average of resistivity among the closest neighbours is computed 
before resistive terms are taken into account in the computation of electromotive 
forces. Resistivity derived in this fashion is then interpolated to the positions of 
the current density components (cell edges) (see Fig. [1} and then the resistive part 
of the electromotive force (£ v = rjj) is computed. Eventually, we compute the 
diffusive terms of induction equation and update magnetic field components. 
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3 Current sheet test 



The following 2D tearing instability setup has been proposed by Hawley & Stone 
(|1995|) (see also the web-page of ATHENA code ( [Stone et al., 2008D ) to test ro- 
bustness of MHD algorithms at extreme conditions. The initial setup assumes 
magnetic field B y = ±1, which changes sign at two locations (x = ±0.25), in a 
domain —0.5 < x < 0.5 and —0.5 < y < 0.5 with doubly periodic boundary con- 
ditions. Initial density and pressure are uniform: p = 1, p = (3/2, where /3 = 0.1 
is the 'plasma-/?'. The initial sinusoidal velocity perturbation is perpendicular to 
magnetic field. The initial setup consists of two current sheets at x = ±0.25, which 
promote magnetic reconnection, due to the numerical or physical resistivity. 

We shall use the above setup to demonstrate qualitative differences between 
simulations with the uniform and current dependent resistivity prescriptions. For 
the uniform resistivity test we assume 771 = 10~ 4 and 772 = 0. It is apparent 
that current sheets form elongated structures, that look rather smooth, since the 
uniform resistivity acts everywhere in the computational domain. The uniform 
resistivity prescription applied in this experiment (Fig. [2]) , leads to configurations 
resembling the Sweet-Parker model (Parker (|1957p ; Sweet (j!958p ). 




Fig. 2. Snapshots of tearing instability at t = 2, 4, 6 with an uniform resistivity. 
Grayscale represents current density and arrows denote normalized magnetic field. Res- 
olution 512 2 

To simulate magnetic reconnection with anomalous resistivity we assume 771 = 0, 
772 = 10 -4 and J cr it = 100. The results of a 2D simulation of the tearing instability, 
with the initial condition specified at the beginning of this section are shown in 
Fig. [3J We find that in the case of localized resistivity multiple magnetic islands 
are formed along each current sheet, and the reconnection process is apparently 
more localized than in the uniform resistivity case. During the evolution smaller is- 
lands undergo coalescence to form bigger ones. Finally, two large magnetic islands 
remain and the reconnection tends to form small-size current sheets, resembling 
'X-type' reconnection points of Petschck's reconnection model ( |Petschek, 1964| ) . 
The reconnection process is faster in the case of localized resistivity, since the final 
volume of magnetic islands is apparently larger than in the uniform resistivity case. 
Although magnetic reconnection can be observed in 'ideal MHD' simulations due 
to numerical resistivity, utilization of the resistivity module allows for controlling 
the speed and properties of the reconnection process in MHD simulations. 
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Fig. 3. Snapshots of tearing instability at t = 0.8, 1.2, 1.6,2,4,6 with anomalous resis- 
tivity. Grayscale, vectors and resolution as in Fig. [2] 
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